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Abstract 

Galam reshuffling introduced in opinion dynamics models is investi- 
gated under the nearest neighbor Ising model on a square lattice using 
Monte Carlo simulations. While the corresponding Galam analytical crit- 
ical temperature Tc ~ 3.09 [J/fcs] is recovered almost exactly, it is proved 
to be different from both values, not reshuffled (Tc = 2/arcsinh(l) « 2.27 
[J/fcs]) and mean-field (Tc = 4 [J/fcs]). On this basis, gradual reshuf- 
fling is studied as function of < p < 1 where p measures the probability 
of spin reshuffling after each Monte Carlo step. The variation of Tc as 
function of p is obtained and exhibits a non-linear behavior. The simplest 
Solomon network realization is noted to reproduce Galam p — 1 result. 
Similarly to the critical temperature, critical exponents are found to differ 
from both, the classical Ising case and the mean field values. 

Keywords: Ising nearest neighbors; reshuffling; sociophysics; phase transition; 
critical temperature. 

1 Introduction 

Sociophysics is based on the use of concepts and tools from physics to describe 
social and political behavior pQ. While the validity of such a transfer has been 
long questioned among physicists 2 , none ever has expected that some basic 
sociophysics question may in turn lead to new development within physics. In 
this paper we report, to our knowledge for the first time, on such a development. 
Indeed economics long ago has contributed substantially to statistical physics, 
but indirectly, providing several tools it developed on its own, which were then 
taken over by physics jjj]. 

Opinion dynamics has been a very active field of research among physicists in 
the recent years |S1 |H] • Most models are based on some local update rule and 
corresponding results are obtained using numerical simulations |SJ . The 
Sznajd model initiated in 2000 [5] is the most popular and has generated a great 
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deal of papers In parallel, a much older one proposed by Galam in 1986 0, 
has stayed confined to only few works [Sj. But a recent universal scheme for 
local rules shows that both models have the same dynamics |S] prompting the 
question of why they have been perceived as different with an advantage for the 
first one. 

Indeed, the reason lies in the fact that, although the Galam model is solvable 
analytically, its very formulation has been perceived as the signature of an in- 
trinsic mean field nature. Since its dynamics is monitored by a total reshuffling 
of agents between repeated local updates, in principle everybody can interact to 
everybody. This fact has been understood as everybody does interact to every- 
body simultaneously, as in a mean field treatment [10) . However that is not the 
case due to the local range of interactions which are restricted to separate small 
groups of agents after each reshuffling. At least, for years, Galam was adamant 
in claiming it. 

In this work we address this question by a thorough numerical Monte Carlo 
investigation of the effect of reshuffling spins on the phase diagram of the two- 
dimensional nearest neighbor (NN) Ising model JT]. Reshuffling is introduced 
gradually according to the variable < p < 1 where p is the probability of 
reshuffling all the spins of the lattice at each Monte Carlo step. We call it the 
Gradually Reshuffled Ising Model and denote it by GRIM. It is worth to stress 
that after each spin reshuffling, interactions stay local among NN. 

The critical temperature Tc is calculated for a series of values of p from 
p = (square lattice Ising model — SLIM) up to p = 1 (totally reshuffled Ising 
model — TRIM). Binder's cumulant for Tc evaluation is used to avoid finite size 
effect ^2] • As expected, at p — the SLIM exact value Tc — 2/arcsinh(l) ss 
2.27[J/fcs] is recovered. The variation of the GRIM Tc as function of p is 
found to exhibit a non-linear behavior. A similar study of gradual reshuffling 
was performed for Galam opinion model in the case of local groups of size four 



The unifying Galam scheme [5] is then used to map the TRIM into a local 
rule model where update groups are of a size five. Its allows an analytical 
calculation of the critical temperature Tc- In contrast to the simulations, it 
is found to depend on the choice of the dynamics. Metropolis and Glauber 
yield different values, respectively Tc ~ 1.59 [J/fcs] and Tc « 3.09[J/fcs]. 
Last one reproduces almost exactly the Monte Carlo result at p — 1, Tc ~ 
3.03[J/fcs]. The simplest realization [Ll of the Solomon network (SN) ^1 is 
noted to reproduce TRIM results. Similarly to the critical temperature, critical 
exponents are found to differ from both, the SLIM case and the mean field 
values. Concluding remarks with respect to future work end the paper. 

2 The gradually reshuffled Ising model 

We consider the nearest neighbor Ising model on a square lattice with 
ferromagnetic interactions 




(i) 



2 



where Si — ±1 is the Ising spin variable at each node i of the square lattice and 



Jij 



J > if i and j are neighbors, 
otherwise. 



is short-range ferromagnetic exchange integral. 

Monte Carlo simulations are performed using either a Glauber or a Metropo- 
lis dynamics. In the Glauber dynamics, every time step all spins at the lattice 
are investigated in type-writer fashion, i.e. spin-by-spin: from left to right and 
then from top to bottom. For each spin i in an initial configuration a new 
configuration rji resulting from the single spin flip Si — > — Si is accepted with a 
probability 

n G = exp(-^/fc B T) 

"h-** exp(-E IH /k B T)+exp(-E Vi /k B T)' [ ' 

where E Vi is the energy of configuration = — E m is the energy of configu- 

ration fii and k B is the Boltzmann constant. When all N spins are investigated 
one Monte Carlo step (MCS) is completed. In Metropolis scheme, the accep- 
tance probability of the new configuration may be expressed in a simple form 

p^ u = min{l,exp[-(^ - E^)/k B T]}. (3) 

But here, at each MCS, all the spins are randomly visited and updated (a 
random list of spins labels assures that each spin is reached exactly once). 

The value of the critical temperature Tc is extracted from the evaluation of 
the order parameter m = Si/N as a function of T. In a second step, to get a 
more precise estimate of Tc in the thermodynamic limit, we calculate Binder's 
fourth cumulant of the order parameter defined by 

where (• • • ) represent the thermal average (taken over 400000 MCS after dis- 
carding 100000 MCS for equilibration at each temperature). This cumulant 
should go to 2/3 below Tc and to zero above Tc when the size increases, and 
the finite-size estimates are expected to cross around Tc 

Our Monte Carlo simulations have been performed on square lattices of size 
50 < L < 500 assuming semi-periodic (helical) or periodic boundary condi- 
tions. For given temperature T, simulation starts with m = 1, i.e., with all 
spins Si — 1. In every MCS, before updating all spins, with a probability p 
the reshuffling procedure takes place: random permutation of all spin labels is 
produced and their positions are rearranged according to that new labels order. 
With probability 1 — p, all spins keep their current position. 



3 The results 

In Fig. m(T) for size L = 500 and different reshuffling probabilities p are 
shown. Also m(T) for SN browsed from Ref. ^1] for a million spins is included. 
The latter agrees very well with reshuffling case for p = 1 (TRIM). In SN a 
random permutation of spin labels is created and each spin interacts with its 
two NN and two additional neighbors of its mirror in reshuffled lattice. In 
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Figure 1: Magnetization m(T) dependence for L — 500 obtained with Glauber 
dynamics. For the SN N = 10 6 spins were simulated. 



contrast to current work the reshuffling procedure takes place only once: before 
the simulation starts and the spin neighborhood is fixed during simulation. The 
spin dynamics is governed via Eqs. (0 and but links are directed, ft 
means, that if spin j is a neighbor of spin i does not yield i being neighbor of 

J- 

Decreasing reshuffling probability p shifts m(T) towards SLIM with Tc — 
2/arcsinh(l) w 2.27 [J/ks] |2Q|- The dependence Tcip) is found to be nonlinear 
as a function of the reshuffling parameter p. In the range p > 0.1 it may be 
approximated with the logarithmic law Tc = 0.19 lnj? + 3.01 (see Fig. |2Jl. 

Examples of the order parameter time evolution m(T) at fixed temperature 
T = 2.5 [J/k B ] are shown in Fig. El for p = (SLIM) and p — 1 (TRIM). The 
curves show that m(T) vanishes for p = as expected from its Tc = 2.27 < 
2.5 [J/fcs]. On the opposite for p = 1 we have Tc = 3.03 > 2.5 [J/fcs] and m(T) 
stabilizes at a non zero value as it should be in the associated ordered phase. 

In Fig. 0|Binder's cumulant U dependence on the temperature T for different 
system sizes L and reshuffling probability p is presented. The common crossing- 
point of these curves predicts Curie temperature Tc- 

In Fig. El critical exponents (3 — which describe magnetization behavior 
in the vicinity of Tc — arc presented. For p = and p — I these exponents 
are respectively Pslim = 0.11 and /3trim = 0.31. The exact value at p = 
is PsliM = V 8 EOIEU while the mean field value is P$ F = 1/2 gTJ. From 
Ptrim = 0.31 ^ 1/2 = /?mfi wc can conclude that TRIM, and thus Galam 
reshuffled models, are not mean field models, since otherwise the associated 
critical exponents would have been mean field exponents, as it must be for any 
mean field approximation. 
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T c =0.19 ln(p)+3.01 




Figure 2: Dependence of Curie temperature Tcip) = 0.19 In p + 3.01 on the 
reshuffling probability p given via cross-point of U(T) curves for different system 
sizes 50 < L < 300. 



N=10 6 , T=2.5 [J/k B ] 



E 




time [MCS] 

Figure 3: Time evolution of the magnetization m(T) obtained with Glauber 
dynamics at T — 2.5 [J/ks] for p — and p = 1. Associated Tc are respectively 
T c = 2.27 [J/k B ] and T c = 3.03 [J/k B ]. Simulations are carried out for N = 10 6 
spins. 
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(a) p=0.5: T c =2.885 [J/k B ] 
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(b) p=1 : T c =3.032 [J/k B ] 
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Figure 4: Binder's cumulant U dependence on the temperature T for different 
system sizes L and reshuffling probability (a) p = 1/2 and (b) p = 1. 
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Figure 5: Critical exponents (3 — which describes magnetization behavior in 
the vicinity of Tc — for p = and p — 1 are /?slim = 0.11 and Ptkim = 0.31, 
respectively. N = 10 6 spins were simulated. Nmcs — 2 • 10 4 and (m) is averaged 
over last 10 4 [MCS]. 

4 The Galam unifying scheme 

We now go one step further in the investigation of the validity of Galam unifying 
scheme 9 . It is a general sequential frame, which operates via local updates of 
small groups of spins randomly selected. Given one group, a specific probabilis- 
tic majority rule is then applied for each peculiar configuration. It is a function 
of the local value of the order parameter, i.e., the ratio of majority to minority. 

Applying this scheme to the square NN Ising model investigated above and 
following the Monte Carlo procedure, which considers one spin and its four NN 
leads to use local groups of one size, with five spins. On this basis, we define <?& 
as the probability that a group of five spins with k plus and (5 — k) minus ends 
up with five plus. Simultaneously the probability to get five minus is (1 — 
With five spins, the number k of plus spins can vary from 5 down to producing 
six different coefficients gk- But from up-down symmetry this number reduces 
to three with go = 1 — g§, g\ = 1 — 54 and g2 — 1 — .93- On this basis, if p(t) is 
the proportion of plus spins at time t, we obtain for the new proportion p(t + 1) 
after one update cycle 



To evaluate the coefficients gt in Eq. © we enumerate all possible 2 5 = 32 
configurations associated to a local group of size five and regroup them by their 
respective number of plus spins. The corresponding energies 2£ Mi are calculated 



^+i) = E(fc)^[^)] fc [ 1 -^)] 5 - fe , 



(•5) 



fc=0 
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as well as energies E Vi are obtained once the central spin has been flipped. As 
described in Sect. 2 the central spin flip is accepted with probability given by 
either Eq. (0) or Q depending on chosen dynamics, Glauber or Metropolis. 
However, here this probability is weighted by the plus spin proportion We 
now illustrate the scheme in one case with four plus spins (Si = +1) and one 
minus (Si = —1). Being on a square lattice, either the minus spin is at the 
center or not. One configuration corresponds to the first case, 



+ 



and for the second case, there exist four configurations of the type 



+ 



where the minus spin can rotate on each one of the four NN sites. 

Using Glauber dynamics, in the first case with one configuration, the flip is 
validated with a probability a weighted by with a proportion one of plus spins. 
In contrast, the initial configuration is conserved with a probability (1 — a) 
weighted by 4/5, the proportion of plus spin. In the second case, the probability 
for a flip is given by 6 weighted with a proportion 3/5 of plus spins while for no 
flip we have respectively (1—6) and 4/5, both values being multiplied by four. 
^From Eq. (0 we have 

exp(AJ/k B T) u _ exp(-2J/fc s T) 



exp(4J/fc s T) + exp(-4J/fc B T) ' exp(2 J/k B T) + cxp(-2 J/k B T) 
Combining above results yields for the coefficient g$ in Eq. (5), 

4 + a (f\ 

In a similar fashion, performing similar above evaluation for the other configu- 
rations we get the coefficients {g{\ with 

20 + a - 46 31-46 
54 = , A go—- (8) 

In parallel, using Metropolis dynamics from Eq. the first configuration 
leads to a flip with probability one (min[l, exp(8J/ k B T)] = 1), and the second 
one is performed with probability 

c = cxp(-4J/fc s T). (9) 

The corresponding coefficients {gi} become 

5-d 21 -4c 28 

55 = — 54 = ^5- and 93 =-, (10) 

where d = exp(— 8J/k B T) is a probability of the configuration change 



+ 
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given by Metropolis dynamics. 

At this stage, we can calculate the unstable fixed points from Eq. J5j to 
get the corresponding critical temperature, respectively Tc = 3.09 [J/ks] for 
Glauber, and Tc — 1.59 [J/fcs] for Metropolis. The Glauber result is rather 
close to the numerical finding Tc = 3.03 [J/ks] but yet a bit different. At 
present we do not understand why. In contrast, the simulation results are in- 
variant in the choice of the dynamics. It is worth to note that a value of Tc a 
little bit larger than 3 [J/ks] was obtained previously by Malarz from a Monte 
Carlo simulation of SN [H] . 

5 Discussion 

In this work, we have studied the effect of gradual spin reshuffling on the phase 
diagram of the lattice square Ising model. Performing a series of Monte Carlo 
simulations, we have proved that Galam reshuffling is neither a mean-field treat- 
ment nor the usual Ising model (SLIM). In particular it yields non mean field 
exponents. 

On this basis, new physical questions arise. First, does reshuffling create a 
new universality class for the Ising model? At which value of the reshuffling 
parameter p does the crossover occur? Besides, the theoretical interest, and the 
sociophysics consequences, to find a physical situation which could correspond 
to reshuffling would be also interesting. 

Acknowledgments. The authors are grateful to Dietrich Stauffer for many 
stimulating discussions and permanent criticism of their ideas, claims and re- 
sults. Part of numerical calculations was carried out in ACK-CYFRONET 
AGH. The machine time on HP Integrity Superdome is financed by the Pol- 
ish Ministry of Science and Information Technology under grant No. MNil/- 
HPJLSD /AGH /047/2004. 



References 



[1] 


D. 


Stauffer, Comput. Sci. Eng. 5, 71 (2003); 




D. 


Stauffer, Fractals 11, 313 (2003); 




D. 


Stauffer, Physica A336, 1 (2004); 




S. 


Galam, Y. Gefen and Y. Shapir, Math. J. Socio. 9, 1 (1982). 


[2] 


S. 


Galam, Physica A336, 49 (2004). 


[3] 


V. 


Pareto, Cours d'Economique Politique, Vol. 2, (Macmillan, Paris, 1897) 


[4] 


D. 


Stauffer, physics/0503115. 




P. 


Gawrohski, P. Gronek, K. Kulakowski, physics/0503085 




P. 


Gawrohski, P. Gronek, K. Kulakowski, physics/0501160 




S. 


Fortunato, cond-mat/0501105 




K. 


Kulakowski, P. Gawrohski, P. Gronek, physics/0501073 



A.O. Sousa, Physica A348, 701 (2005); 



9 



S. Fortunato, Physica A348, 683 (2005); 
S. Fortunato, cond-mat/0405083 

E. Ben-Nairn, Europhys. Lett. 69, 671 (2005); 

M. F. Laguna, S. Risau-Gusman, G. Abramson, S. Goncalves, J. R. Iglesias, 
nlin. AO/0404024; 

M.F. He, Q. Sun, H.S. Wang, Int. J. Mod. Phys. C15, 767 (2004); 
D. Stauffer, A.O. Sousa, C. Schulze, cond-mat/0310243 

C. J. Tessone, R. Toral, P. Amengual, H.S. Wio, M. San Miguel, Eur. Phys. 
J. B39, 535 (2004); 

D. Stauffer, H. Meyer-Ortmanns, Int. J. Mod. Phys. C15, 241 (2004); 

R. Hegselmann, U. Krause, J. Artificial Soc. Soc. Simul. 4 5 (3) 2 (2002); 

G. Deffuant, D. Neau, F. Amblard, G. Weisbuch, Complex Syst. 3, 87 
(2000); 

G. Weisbuch, G. Deffuant, F. Amblard, J.-P. Nadal, Complexity 7, 55 
(2002); 

G. Deffuant, F. Amblard, G. Weisbuch, T. Faure, J. Artificial Soc. Soc. 
Simul. 5 (4) 1 (2002); 

J. A. Holyst, K. Kacperski, F. Schweitzer, Ann. Rev. Comput. Phys. IX, 
253 (2001) and references therein. 

[5] K. Sznajd-Weron, J. Sznajd, Int. J. Mod. Phys. Cll, 1157 (2000). 

[6] C. Schulze, Int. J. Mod. Phys. C15, 867 (2004); 

M.F. He, B. Li, L.D. Luo, Int. J. Mod. Phys. C15, 997 (2004); 

D. Stauffer, J. Math. Sociol. 28, 25 (2004); 

L. Sabatclli, P. Richmond, Physica A334, 274 (2004); 

J.J. Schneider, Int. J. Mod. Phys. C15, 659 (2004); 

C. Schulze, Int. J. Mod. Phys. C15, 569 (2004); 

L. Behera, F. Schweitzer, Int. J. Mod. Phys. C14, 1331 (2003); 

L. Sabatelli, P. Richmond, Int. J. Mod. Phys. C14, 1223 (2003); 

J. Bonnekoh, Int. J. Mod. Phys. C14, 1231 (2003); 

F. Slanina, H. Lavicka, Eur. Phys. J. B35, 279 (2003); 
C. Schulze, Physica A324, 717 (2003); 

C. Schulze, Int. J. Mod. Phys. C14, 95 (2003); 

D. Stauffer, P.M.C. de Oliveira, Eur. Phys. J. B30, 587 (2002); 
D. Stauffer, Comput. Phys. Commun. 146, 93 (2002); 

D. Stauffer, Int. J. Mod. Phys. C13, 315 (2002); 

A.T. Bernardes, D. Stauffer, J. Kertesz, Eur. Phys. J. B25, 123 (2002); 
A.S. Elgazzar, Int. J. Mod. Phys. C12, 1537 (2001); 

4 available online at http://jasss.soc.surrey.ac.uk/ 



10 



I. Chang, Int. J. Mod. Phys. C12, 1509 (2001); 

R. Ochrombel, Int. J. Mod. Phys. C12, 1091 (2001); 

A.T. Bernardes, U.M.S. Costa, A.D. Araujo, D. Stauffer, Int. J. Mod. Phys. 
C12, 159 (2001); 

A. A. Moreira, J.S. Andrade, D. Stauffer, Int. J. Mod. Phys. C12, 39 (2001); 

D. Stauffer, A.O. Sousa, S.M. de Oliveira, Int. J. Mod. Phys. Cll, 1239 
(2000). 

[7] S. Galam, J. Math. Psychology 30, 426 (1986); 
S. Galam, J. Stat. Phys. 61, 943 (1990); 
S. Galam, Physica A285, 66 (2000). 

[8] S. Galam, Eur. Phys. J. B25, 403 (2002). 

[9] S. Galam, [cond-mat/ 0409484 

[10] C.J. Tessone, R. Toral, P. Amengual, H.S. Wio, M. San Miguel, Eur. Phys. 
J. B39, 535 (2004). 

[11] W. Lenz, Z. Phys. 21 (1920) 613; 

E. Ising, Z. Phys. 31 (1925) 253. 

[12] K. Binder, D.W. Heermann, Monte Carlo Simulation in Statistical Phyics, 
(Springer Verlag, 1988). 

[13] S. Galam, B. Chopard, A. Masselot, M. Droz, Eur. Phys. J. B4, 529 (1998). 

[14] K. Malarz, Int. J. Mod. Phys. C14, 561 (2003). 

[15] T. Erez, M. Hohnisch, S. Solomon, Economics: Complex Windows, Eds. 
M. Salzano and A. Kirman, (Springer, 2005), p. 201, |cond-mat7 0406369 . 

[16] M.C. Gonzalez, A.O. Sousa, H.J. Herrmann, Int. J. Mod. Phys. C15, 45 
(2004). 

[17] D. Stauffer, H. Meyer-Ortmanns, Int. J. Mod. Phys. C15, 241 (2004). 
[18] F. Slanina, H. Lavicka, Eur. Phys. J. B35, 279 (2003). 
[19] F. Schweitzer, J. Zimmermann, H. Muhlenbcin, Physica A303, 189 (2002). 
[20] L. Onsager, Phys. Rev. 65, 117 (1944). 

[21] J.J. Binney, N.J. Dowrick, A.J. Fisher, M.E.J. Newman, A theory of criti- 
cal phenomena. An Introduction to the renormalization group. (Clarendon 
Press, Oxford, 1992) 



11 



